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Abstract 

(N . 

Taking advantage of the recent literature on exact simulation algorithms (Beskos et al. [I]) and 
unbiased estimation of the expectation of certain functional integrals (Wagner 23 , Beskos et al. [5] and 
Fearnhead et al. [6]), we apply an exact simulation based technique for pricing continuous arithmetic 
average Asian options in the Black & Scholes framework. Unlike existing Monte Carlo methods, we are 
no longer prone to the discretization bias resulting from the approximation of continuous time processes 
through discrete sampling. Numerical results of simulation studies are presented and variance reduction 
problems are considered. 
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Introduction 

m 

Although the Black & Scholes framework is very simple, it is still a challenging task to efficiently price 
| Asian options. Since we do not know explicitly the distribution of the arithmetic sum of log-normal variables, 

. there is no closed form solution for the price of an Asian option. By the early nineties, many researchers 

attempted to address this problem and hence different approaches were studied including analytic approxima- 
tions (see Turnball and Wakeman [3U], Vorst [22, Levy [TS] and more recently Lord [IB]), PDE methods (see 
Vecer [2T] . Rogers and Shi 18 , Ingersoll Dubois and Lelievre [5J), Laplace transform inversion methods 
(see Geman and Yor [TO], Geman and Eydeland [5]) and, of course, Monte Carlo simulation methods (see 
Kemna and Vorst [13], Broadie and Glasserman [3], Fu et al. [7]). 

Monte Carlo simulation can be computationally expensive because of the usual statistical error. Variance 
reduction techniques are then essential to accelerate the convergence (one of the most efficient techniques is 
the Kemnafe Vorst control variate based on the geometric average). One must also account for the inherent 
discretization bias resulting from approximating the continuous average of the stock price with a discrete 
one. It is crucial to choose with care the discretization scheme in order to have an accurate solution (see 
Lapeyre and Temam [2]). The main contribution of our work is to fully address this last feature by the use, 
after a suitable change of variables, of an exact simulation method inspired from the recent work of Beskos 
et al. [U [5] and Fearnhead et al. [BJ . 

In the first part of the paper, we recall the algorithm introduced by Beskos et al. [I] in order to simulate 
sample-paths of processes solving one-dimensional stochastic differential equations. By a suitable change of 
variables, one may suppose that the diffusion coefficient is equal to one. Then, according to the Girsanov 
theorem, one may deal with the drift coefficient by introducing an exponential martingale weight. Because 
of the one-dimensional setting, the stochastic integral in this exponential weight is equal to a standard 
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integral with respect to the time variable up to the addition of a function of the terminal value of the 
path. Under suitable assumptions, conditionally on a Brownian path, an event with probability equal to 
the normalized exponential weight can be simulated using a Poisson point process. This allows to accept or 
reject this Brownian path as a path solution to the SDE with diffusion coefficient equal to one. In finance, 
one is interested in computing expectations rather than exact simulation of the paths. In this perspective, 
computation of the exponential importance sampling weight is enough. The entire series expansion of the 
exponential function permits to replace this exponential weight by a computable weight with the same 
conditional expectation given the Brownian path. This idea was first introduced by Wagner [231 (Ml [55] 
in a statistical physics context and it was very recently revisited by Beskos et al. [2] and Fearnhead et al. 
[B] for the estimation of partially observed diffusions. Some of the assumptions necessary to implement the 
exact algorithm of Beskos et al. [I] can then be weakened. 

The second part is devoted to the application of these methods to option pricing within the 

Black & Scholes framework. Throughout the paper, St — So exp (aWt + (r — 8 — "^-)* J represents the 

stock price at time t, T the maturity of the option, r the short interest rate, a the volatility parameter, S 
the dividend rate and (M^)tg[o,T] denotes a standard Brownian motion on the risk- neutral probability space 

(fi,J-", P). We are interested in computing the price Co = E (e~ rT f (cxSt + P Jq Stdtj^j of a European 

option with pay-off / (cxSt + fi Jq StdtJ assumed to be square integrable under the risk neutral measure P. 
The constants a and j3 are two given non-negative parameters. 

When a > 0, we remark that, by a change of variables inspired by Rogers and Shi [TS], aSx + Jq Stdt 
has the same law as the solution at time T of a well-chosen one-dimensional stochastic differential equation. 
Then it is easy to implement the exact methods previously presented. The case a = of standard Asian 
options is more intricate. The previous approach does not work and we propose a new change of variables 
which is singular at initial time. It is not possible to implement neither the exact simulation algorithm nor 
the method based on the unbiased estimator of Wagner |23j and we propose a pseudo-exact hybrid method 
which appears as an extension of the exact simulation algorithm. In both cases, one first replaces the integral 
with respect to the time variable in the function / by an integral with respect to time in the exponential 
function. Because of the nice properties of this last function, exact computation is possible. 



1 Exact Simulation techniques 

1.1 The exact simulation method of Beskos et al. [1] 

In a recent paper, Beskos et al. I proposed an algorithm which allows to simulate exactly the solution 
of a 1-dimensional stochastic differential equation. Under some hypotheses, they manage to implement an 
acceptance-rejection algorithm over the whole path of the solution, based on recursive simulation of a biased 
Brownian motion. Let us briefly recall their methodology. We refer to [T] for the demonstrations and a 
detailed presentation. 

Consider the stochastic process (£t)o<t<T determined as the solution of a general stochastic differential 
equation of the form : 

f d& = b(£ t )dt + <T(Z t )dW t m 

where b and a are scalar functions satisfying the usual Lipschitz and growth conditions with a non vanishing. 
To simplify this equation, Beskos et al. 1 suggest to use the following change of variables : Xt = rj(£t) 
where r\ is a primitive of ^ (r)(x) = f x -^nj^du). 



2 



Under the additional assumption that — is continuously differentiable, one can apply Ito's lemma to get 



dXt 



v'&Ht + -r,"(Zt)d<Z,Z> t 
b -^dt + dW t -^dt 



dt + dW t 



a(X t ) 



So £t = 77 1 (Xt) where (Xt)t is a solution of the stochastic differential equation 

dX t = a{X t )dt + dW t 
X = x. 



(2) 



Thus, without loss of generality, one can start from equation ^ instead of |T]). 

Let us denote by (Wt)tt[o,T\ the process (W t +^)t£[o,T]i by its law and by Qx the law of the process 
(X t )t£[o.T]- From now on, we will denote by (^t)te[o,r] the canonical process, that is the coordinate mapping 
on the set C([0,T], R) of real continuous maps on [0,T] (see Revuz and Yor [17] or Karatzas and Shreve 

One needs the following assumption to be true 
Assumption 1 : Under Qvy x 1 the process 



L t = exp 



/ a(Y u )dY u - ~ / a 2 (Y u )dv 
Jo z Jo 



is a martingale. 



According to Rydberg [T5] (see the proof of Proposition 2] where we give his argument on a specific 
example), a sufficient condition for this assumption to hold is 
-Existence and uniqueness in law of a solution to the SDE <j2j) - 



-V* E [0,T], f a 2 (Y u )du < 00, Q x and Q w * almost surely on C([0,T], R). 
Jo 



Thanks to this assumption, one can apply the Girsanov theorem to get that Qx is absolutely continuous 
with respect to and its Radon-Nikodym derivative is equal to 



d$x 



exp 



a{Y t )dY t 



1 



a 2 (Y t )dt 



Consider A the primitive of the drift a, and assume that 
Assumption 2 : a is continuously differentiable. 

Since, by Ito's lemma, A(W|) = A(a:) + J Q T ^Wf^W^ + § J Q T a'(W?)dt, we have 



r/Q 



dQw* 



exp 



A(T T ) - A(ar) - ~ / a 2 {Y t ) + a'(Y t )dt 
1 Jo 



Before setting up an acceptance-rejection algorithm using this Radon-Nikodym derivative, a last step is 

needed. To ensure the existence of a density h{u) proportional to exp(A(u) — 2 j ), it is necessary and 
sufficient that the following assumption holds 
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Assumption 3 : The function u H ► exp(A(w) — ^^-) is integrable. 

Finally, let us define a process Z t distributed according to the following law Qz 

Qz = j £((Wr) te [o,r]|W| = y)h(y)dy 

where the notation £(.|.) stands for the conditional law. One has 



d€lz dQ w * dQz 



C exp 



a 2 (Y t ) + a'(Y t )dt 



where C is a normalizing constant. At this level, Beskos et al. [T] need another assumption 



Assumption 4 '■ The function (j) : x t— > 



a 2 (x)-\-a f (x) 



is bounded from below. 



Therefore, one can find a lower bound k of this function and eventually the Radon-Nikodym derivative 
of the change of measure between X and Z takes the form 



dQx 
d^z 



Ce 



-kT 



exp 



cj)(Y t ) -kdt 



The idea behind the exact algorithm is the following : suppose that one is able to simulate a continuous 
path Zt{uj) distributed according to Qz and let M(lo) be an upper bound of the mapping 1 1— > (f>(Z t (uj)) — k. 
Let N be an independent random variable which follows the Poisson distribution with parameter TM(lo) and 
let (Ui, Vi)i-i,,,N be a sequence of independent random variables uniformly distributed on [0, T] x [0, M{uj)\. 
Then, the number of points (Ui, Vi) which fall below the graph {(t, (f>(Z t (ui)) — k);t £ [0, T]} is equal to zero 
with probability exp — J^(f>(Z t (uj)) — kdt . Actually, simulating the whole path [Zt)te[a.T] is not necessary. 

It is sufficient to determine an upper bound for <j>(Z t ) — k since, as pointed out by the authors, it is possible 
to simulate recursively a Brownian motion on a bounded time interval by first simulating its endpoint, then 
simulating its minimum or its maximum and finally simulating the other pointt[^|. For this reason, one needs 
the following assumption for the algorithm to be feasible : 

Assumption 5 : Either lim sup 4>(u) < +oo or limsup0(u) < +oo. 

u— >+oo u— >- — oo 

Suppose for example that limsup0(u) < +oo. The exact algorithm of Bekos et al. 1 then takes the 

U— ¥ + 00 

following form : 
Algorithm 1 

1. Draw the ending point Zt of the process Z with respect to the density h. 

2. Simulate the minimum m of the process Z given Zt ■ 

3. Fix an upper bound M (to) = sup{</)(it) — fc; u > to} for the mapping 1 1— >• 4>(Z t ) — k. 

4- Draw N according to the Poisson distribution with parameter TM(m) and draw (Ui, T^)i=i...jv> o. 
sequence of independent variables uniformly distributed on [0,T] x [0, M(m)]. 

5. Fill in the path of Z at the remaining times (Ui)i = \...N ■ 



J In their paper, the authors explain how to do such a decomposition of the Brownian path. 
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6. Evaluate the number of points (^i)i=i...Ar such that Vi < <$>(Zu^) — 
If it is equal to zero, then return the simulated path Z . 
Else, return to step 1. 



k. 



This algorithm gives exact skeletons of the process X, solution of the SDE ©. Once accepted, a path can 
be further recursively simulated at additional times without any other acceptance/rejection criteria. We also 
point out that the same technique can be generalized by replacing the Brownian motion in the law of the 
proposal Z by any process that one is able to simulate recursively by first simulating its ending point, its 
minimum/maximum and then the other points. Also, the extension of the algorithm to the inhomogeneous 
case, where the drift coefficient a in ([5]), and therefore the function <j), depend on the time variable t, is 
straightforward given that the assumptions presented above are appropriately modified. 



1.2 The unbiased estimator (U.E) 

In finance, the pricing of contingent claims often comes down to the problem of computing an expectation 
of the form 

C = E(/(Xr)) (3) 

where X is a solution of the SDE ([2]) and / is a scalar function such that /(-XV) is square integrable. In a 
simulation based approach, one is usually unable to exhibit an explicit solution of this SDE and will therefore 
resort to numerical discretization schemes, such as the Euler or Milstein schemes, which introduce a bias. 
Of course, the exact algorithm presented above avoids this bias. Here, we are going to present a technique 
which permits to compute exactly the expectation ([3]) while assumptions 4 and 5 on the function a which 
appears in the Radon-Nikodym derivative are relaxed. 

Using the previous results and notations, we get, under the assumptions 1 and 2, that 



C = E /(W^)exp 



A(W%) - A(x) 



1 



a 2 (W t x ) + a'(W t x )dt 



(4) 



In order to implement an importance sampling method, let us introduce a positive density p on the real 
line and a process (Z t )t^[o,T] distributed according to the following law Qz 



By one has 



C = EU(2 T )exp 



4>{z t )dt 



(5) 



where ip : z t— > f(z)' 



, ; 2T — and 6 : z n- £LJ£i±£_i£2, We do not impose p to be equal to the 
density h of the previous section. It is a free parameter chosen in such a way that it reduces the variance of 
the simulation. 

In his first paper, Wagner |23j constructs an unbiased estimator of the expectation ([5]) when i/j is a 
constant, (Zt)te[o,T] is an valued Markov process with known transition function and tj> is a measurable 

function such that E < +oo. His main idea is to expand the exponential term in a power 

series, then, using the transition function of the underlying Markov process and symmetry arguments, he 
constructs a signed measure v on the space y — \Jn= o ([0,T] x R d ) n+1 such that the expectation at hand is 
equal to v(y). Consequently, any probability measure /ionY that is absolutely continuous with respect to 
v gives rise to an unbiased estimator Q defined on (y,fi) via C(y) = j^(y)- In practice, a suitable way to 
construct such an estimator is to use a Markov chain with an absorbing state. Wagner also discusses variance 
reduction techniques, specially importance sampling and a shift procedure consisting on adding a constant c 
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to the integrand <f> and then multiplying by the factor e~ cT in order to get the right expectation. Wagner [25] 
extends the class of unbiased estimators by perturbing the integrand <fi by a suitably chosen function </>o and 
then using mixed integration formulas representation. Very recently, Beskos et al. [5] obtained a simplified 
unbiased estimator for ([5]), termed Poisson estimator, using Wagner's idea of expanding the exponential in 
a power series and his shift procedure. To be specific, the Poisson estimator writes 

^Z T )e^f[ C -^l (6 ) 

1=1 

where N is a Poisson random variable with parameter cp and (Vj), is a sequence of independent random 
variables uniformly distributed on [0,T]. Fearnhead et al. [6] generalized this estimator allowing c and cp to 
depend on Z and N to be distributed according to any positive probability distribution on N. They termed 
the new estimator the generalized Poisson estimator. We introduce a new degree of freedom by allowing the 
sequence (Vi)i to be distributed according to any positive density on [0,T]. This gives rise to the following 
unbiased estimator for (J5| : 



Lemma 1 — Letpz and qz denote respectively a positive probability measure on IKI and a positive probability 
density on [0,T]. Let N be distributed according to pz and (Vi)j£iM» be a sequence of independent random 
variables identically distributed according to the density qz, both independent from each other conditionally 
on the process (Z t )te[o,T]- Let cz be a real number which may depend on Z. Assume that 



E I \^(Z T )\e- CzT exp 



\c z - <t>{Z t )\dt 



< oo. 



The 



Pz(N)Nl£± qz{Vi) 



is an unbiased estimator of Co . 



Proof : The result follows from 



E MZ T )e 



-c z T 



N 

n 



c z - (j>{Z Vi ) 



p z (7V)7V!ii q z {Vi) 



\ T ^(loCz-HZt)dt) n 
(Z t ) t ^ T] )^^(Z T )e^ T J2 S 



Pz{n) n! 



Pz(n) 



*p(Z T )exp - / <f>{Z t )dt 



Using ([7]). one is now able to compute the expectation at hand by a simple Monte Carlo simulation. The 
practical choice of pz and qz conditionally on Z is studied in the appendix 14.11 

As pointed out in Fearnhead et al. [6], this method is an extension of the exact algorithm method since, 
under assumptions 3, 4 and 5, the reinforced integrability assumption of Lemma [T] is always satisfied. 

Indeed, suppose for example that lim sup (/>(it) < +oo and let k be a lower bound of <j>, mz be the 

minimum of the process Z and M z an upper bound of {4>{u) — k,u > mz}. Then, taking cz — Mz + k in 



(i 



Lemma [T] ensures the integrability condition : 

E ^j(Z T )\e- ( - Mz+ ^ T eIo lMz+k -' f ' {Zt)ldt ^ = E ^(Z T )\er^ Mz+k ^ T ei» Mz+k ~^ z ^ dt 

= E(\ip(Z T )\e-f° ^ dt> ) <oo 

and hence, one is allowed to write that 

Co = E f ^ T )e-(^+^^I— n Mz + k .^A . 
y pz{N)Nl gz(n) y 

Better still, the random variable -0(Zr)e~^ Mz+fc - )T pz ^ N \ IliLi ~^iT^7~^ * s sc l uare integrable when 
is the Poisson distribution with parameter MzT + A: and is the uniform distribution on [0, T] since we 
have then 



< 



E(^ 2 (^t)) < 00. 



The last inequality follows from the square integrability of / : whenever one is able to simulate from the 
density h, introduced in the exact algorithm, by doing rejection sampling, there exists a density p such that 
ip, which is equal to f(Zx) ^[zt) u P ^° a cons t an t factor, is dominated by / and so is square integrable. 

The square integrability property is very important in that we use a Monte Carlo method. We see that, 
whenever the exact algorithm is feasible, the unbiased estimator of lemma [T] is a simulable square integrable 
random variable, at least for the previous choice of pz and qz- 

Remark 2 — One can derive two estimators of Co from the result of Lemma [7] : 

1" . e A(Z^-A( x) -^l ^ t J N* C Z - ftZty 

^ - n^ {Z%T) V^pJW) 6 " PzmmU qz (VJ) 3 

n ^ eA{ z' T )-A( x )-^l 1 ^Cz-t(Zy 



2^ P (z* T ) p z (N i )m\f = \ qz {y;) 
{-[ V^pjzf) p~zWW\ 1 = 1 q z iyy) 



2 Application : the pricing of continuous Asian options 

In the Black & Scholes model, the stock price is the solution of the following SDE under the risk-neutral 
measure P 

^ = (r - S)alt + adW t (8) 

where all the parameters are constant : r is the short interest rate, S is the dividend rate and a is the 
volatility. 

Throughout, we denote 7 = r — 8 — The path- wise unique solution of © is 

S t = S exp(o-W t + -ft) . 
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We consider an option with pay-off of the form 

/ LtSr + 13 £ S t dtj (9) 

where / is a given function such that E (j 2 (aSr + P Jq S t dtj^j < oo, T is the maturity of the option and 

a, /3 are two given non negative parameter^. Note that for a = 0, this is the pay-off of a standard continuous 
Asian option. 

The fundamental theorem of arbitrage- free pricing ensures that the price of the option under consideration 

is 

C = E ( e- rT f \aS T + P I S u du\ ) . 



At first sight, the problem seems to involve two variables : the stock price and the integral of the stock 
price with respect to time. Dealing with the PDE associated with Asian option pricing, Rogers and Rogers 
and Shi [TS] used a suitable change of variables to reduce the spatial dimension of the problem to one. We 
are going to use a similar idea. 

Let 



We have that 



aS e aBt+ ^ +/3S* f e aB ^ s ds 
Jo 



where we set B s — W t — W t ~ s , Vs £ [0, i]. Clearly, (B s ) s& i 0tt ] is a Brownian motion and thus the following 
lemma holds 

r* 

Lemma 3 — Vf G [0, T], £ 4 and aSt + P \ S u du have the same law. 

Jo 

As a consequence 

Co = E {e- rT f{H T )) . 

By applying Ito's lemma, we verify that the process (£t)t>o is a positive solution of the following 1-dimensional 
stochastic differential equation for which path-wise uniqueness holds 

| d£t = pS dt + Z t (adWt + (~/+^)dt) 
{ £o = aSo. 

We are thus able to value Co by Monte Carlo simulation without resorting to discretization schemes using 
one of the exact simulation techniques described in the previous section. In the case a = 0, one has to deal 
with the fact that £ t starts from zero which is the reason why we distinguish two cases. 



3 The underlying of this option is a weighted average of the stock price at maturity and the running average of the stock 
price until maturity with respective weights a and 0T. 
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2.1 The case a ^ 

We are going to apply both the exact algorithm of Beskos et al. 1 and the method based on the unbiased 
estimator of lemma [T] 

We make the following change of variables to have a diffusion coefficient equal to 1 : 

° 1 X = x with x= iSS^o). 

Thus 

C = E(e- rT f(e" x T)). 
The following proposition ensures that assumption 1 is satisfied. 



Proposition 4 — The process (£{)t£[o,T] defined by 
L t = exp 

is a martingale under Qw* 



Q a a 2 J a a 



Proof : Under Qw, (Lt)te[o,T\ is clearly a non-negative local martingale and hence a super-martingale. 
Then, it is a true martingale if and only if Eq w . x (Lt) = 1- 

Checking the classical Novikov's or Kamazaki's criteria is not straightforward. Instead, we are going to 
use the approach developed by Rydberg [19] (see also Wong and Heyde [27]) who takes advantage of the link 
between explosions of SDEs and the martingale property of stochastic exponentials. 

Let us define the following stopping times : 

T n (Y) = inf It e R + such that f (- + ^-e'^A du > n\ , 



o 



with the convention inf{0} = +oo. 

The stopped process (L tATn ry))t&[o,T] is a true martingale under Qw* since Novikov's condition is fulfilled. 
According to the Girsanov theorem, one can define a new probability measure Q^-, which is absolutely 
continuous with respect to Qw* > by its Radon-Nikodym derivative 



x 



dQ\yic 



L 



TAr n {Y)- 



Hence 

EQJ ( 1 {r„(Y)>T}) = Eq w x (l{r„(y)>T}L T Ar„(y)) • 

Since (T n (Y))„^ is a non decreasing sequence, we can pass to the limit in the right hand side We get 
where ^(Y) denotes the limit of the non decreasing sequence (T n (Y)) ne ^. 

Under Qw*, (^*)te[o,T] has the same law as a Brownian motion starting from x so t 00 (Y) = +oo , Qw x 
almost surely, and consequently 

Eq w ,(Lt)= hm €l n x (r n (Y) > T) . 
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On the other hand, the Girsanov theorem implies that, under Q^-, (Yt)te[o,TAr„(Y)] solves a SDE of the 
form (fTTj) . To conclude the proof, it is sufficient to check that trajectorial uniqueness holds for this SDE. 
Indeed, the law of {Yt)te[o,TAT„(Y)] under Q x is the same as the law of (^ / t)tG[o,TAr„(y)] under Qx- Hence 



£&(r n (Y) >T) = Qx(r„0n > T) — > Q x { To0 {Y ) > T) . 

n— f+oo 



Clearly, /„' ( 



2 _i_ ggpg-o-n 



du < co, Qx almost surely, so 

E q ^(L T )=Qx (r 00 (Y)>T) = l 



as required. 

In order to check trajectorial uniqueness for the SDE (fTTj) , we consider two solutions X 1 and X 2 . Wc 
have that 



d{Xl - X?) = ^ (e-^ 1 - e -*?) dt => d\X] - X?\ = ^sign{Xt - X?) ( e -*i - 



dt. 



So 



\X} X}\ = &°f sign{Xl - Xf) ( e -** 
° Jo v 



- e - aX n ds < o. 



The last inequality follows from the fact that x i— > e 51 is a decreasing function. Finally, almost surely, 
Vi > 0, X\ = X? which leads to strong uniqueness. □ 



Consequently, thanks to the Girsanov theorem, we have 



r/Q 



x 



dQw x 



exp 



Jo „ ° ° , 2 Jo CT (T 



Set = J n a(a;)da; = > + ^(1 - e~ au ). Then 



x 



dQw=° 



exp 



1 



A(y r ) - A(x) - ~ / a 2 (F 4 ) + a'(F 4 )^ 
* Jo 



The function u i-> exp (a(u) - (n ^ o) ^ = exp + ^(1 - e"*™) 
can define a new process (Z t ) t6 [ .T] distributed according to the following law Qz 



(12) 



Y °- > ' is clearly integrable so we 



Qz= / c((Wt) te[ o,T]\WT = y)h{y)dy 

J [R 

where the probability density h is of the form 

h(u) = C exp ^A(ii) — - — ) wn ^ ^ a norrnanzm g constant. 



(13) 



Remark 5 — Simulating from this probability distribution is not difficult (see the appendix \ j.2\ for an 
appropriate method of acceptance/rejection sampling). 
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We have 



Cexp 



-(a 2 (Y t )+a'(Y t ))dt 



Set fa) = = G + ^e-™)'-es e- 



inf <b{x) 
xeR 



Set k = inf^giR 4>(x). Finally, we get 



A direct calculation gives 
if 2 7 > a 2 
^(i lo g(^)) otherwise. 



2 " 

2rr 



Ce- fcT exp 



</>(Y t ) - Jfedt 



We check that 



lim 4>(x) = — — < oo 

x— s-+oo 2(7 

lim = +00. 



Hence we can apply the algorithm [T] to simulate exactly Xt and compute Co = E (e rT f(e aXT )) by 
Monte Carlo. On the other hand, using (|12p we get 



1 



- - - / a 2 (W?) + a'(W?)dt 



Cv = E\^- rT f(e° w T)exp 

and we can also use the unbiased estimator presented in the previous section to compute this expectation. 



Remark 6 — We also applied the exact algorithm based on a geometric Brownian motion instead of the 
standard Brownian motion which seems more intuitive given the form of the SDE (Q7J). The algorithm is 
feasible because we can simulate recursively a drifted Brownian motion and therefore a geometric Brownian 
motion by an exponential change of variables. The results we obtained were not different from the first 
method. 



2.1.1 Numerical computation 

For numerical tests, we consider the case 

f(x) = (x-K)+ 

which corresponds to the European call option with strike K . Using the exact simulation algorithm presented 
above, we can simulate the underlying aSx + P Jq Stdt at maturity (see Figured]). Then, all we have to 
do is a simple Monte Carlo method to get the price of the option under consideration. Using the unbiased 
estimator, we get 



^MW-IS^ , Mz+k)T 1 " M Z + k - 4>{Z Vi ) 



Co = E e- rT (e° z r - K)+- = — e^ M ^ T * , TT 

\ V2^p(Z T ) Pz (N)N\}-_\ 



2tt P (Zt) Pz(N)N\ 11 q(Vi) 

where (Z t )te[o,T], Mz, k,pz and qz are defined as in section [PI In order to ensure square integrability, 
we choose pz to be a Poisson distribution with parameter MzT + k and qz to be the uniform distribution on 
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[0, T]. For the density p, a good choice is to consider the density that we use to simulate from the distribution 
h by rejection sampling. 

We test these exact methods against a standard discretization scheme with the variance reduction tech- 
nique of Kemna and Vorst [12] . As pointed out by Lapeyre and Temam [T?] , the discretization of the integral 
by a simple Riemannian sum is not efficient. Instead, we use the trapezoidal discretization. In the sequel, we 
will denote this method by Trap+KV. The table [TJ gives the results we obtained for the following arbitrary 
set of parameters : S = 100, K = 100, r = 0.05, a = 0.3, 5 = 0, T = 1, a = 0.6 and = 0.4. The 
computation has been made on a computer with a 2.8 Ghz Intel Penthium 4 processor. We intentionally 
choose a large number of simulations in order to show the influence of the number of time steps when using 
a discretization scheme. 



Exact Simulation of the underlying : cxSt + /3 / Stdt 



0.016 



0.014 



0.012 




0.000 



St pT 

cxSt + ft I Stdt 
Jo 



270 



313 



357 



400 



Figure 1: Histogram of 10 5 independent realizations of o.St + P L Stdt for a = 0.6 and j3 — 0.4 compared 
with the lognormal distribution of St- 



Method 


M 


N 


Acceptance rate 


Price 


C.I at 95% 


CPU 


Trap+KV 


10 


10 6 




11.46 


[11.43,11.48] 


5 s 


20 


11.46 


[11.43,11.49] 


9 s 


50 


11.47 


[11.44,11.5] 


21 s 


Exact Simulation 




lO 6 


24% 


11.46 


[11.43,11.5] 


81 s 


U.E (c P = M z ,c z = Mz + k) 




10 6 




11.46 


[11.43,11.49] 


17 s 


U.E (c P =c z = l/T) 




lO 6 




11.46 


[11.43,11.49] 


6 s 



Table 1: Price of the option (|9]) using a standard discretization technique and exact simulation methods. 
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Empirical evidence shows that the exact simulation method is quite slow. This is mainly due to the 
fact that the rejection algorithm has a little acceptance rate (24% according to table [T]). Using a geometric 
Brownian motion instead of a standard Brownian motion did not improve the results. Also, simulating 
recursively a Brownian path conditionally on its terminal value and its minimum is time consuming. 

The unbiased estimator is more efficient, especially when we can avoid the recursive simulation of the 
Brownian path. To do so, we choose for pz a Poisson distribution with mean cpT where cp is a free 
parameter. If we assume that the integrability condition in lemma [T] holds, then we can write that 



Co = E I e~ rT {e a '^ - K) 



a A(Z T )-A(x)- 



a c P T-c z T 



2np(Z T ) 



N 

n 

i=l 



CZ 



Regarding the dependence of the exact simulation method with respect to the parameters a and /3, it 
is intuitive that whenever a » /3, the method performs well since the logarithm of the underlying is not 
far from the logarithm of the geometric Brownian motion on which we do rejection-sampling. The table [5] 
confirms this intuition. We see that we cannot apply the algorithm for small values of a and then let a — > 
to treat the case a = 0. 



a + f3 


0.3 


0.4 


0.5 


0.6 


0.7 


Acceptance Rate 


0.003% 


0.47% 


5.66% 


24.43% 


53.85% 



Table 2: Influence of the parameter -^-a on the acceptance rate of the exact algorithm. 



2.2 Standard Asian options : the case a = and (3 > 

A standard Asian option is a European option on the average of the stock price over a determined period 
until maturity. An Asian call, for example, has a pay-off of the form (y S u du — K) + - With our previous 
notations, it corresponds to the case a = 0, f3 — ^ and f(x) = (x — K) + . 

The change of variables we used above is no longer suitable because it starts from zero when a = 0. 
Instead, we consider the following new definition of the process £ 



1 & = So. ° 



(14) 



Obviously, the two variables £t and y f Q S u du have the same law. Hence, the price of the Asian option 
becomes 



Co = E e~ rT 



f^£s u duj \ =E(e- rT f($ T )) 



Remark 7 — The pricing of floating strike Asian options is also straightforward using this method. It is 
even more natural to consider these options since it unveils the appropriate change of variables as we shall 
see below. 

Let us consider a floating strike Asian call for example. We have to compute 



i J S u du - S T )^j ■ 



13 



Using St = Ste as a numeraire (see the seminal paper of Geman et al. 19]), we immediately obtain that 

C a = E Ps Ue^(ir^ d u-l) + ) 



where Pg is the probability measure associated to the numeraire St- It is defined by its Radon-Nikodym 
derivative = e °' 1/lA r-V' r . 

Under Pjr, the process B t = W t — at is a Brownian motion and we can write that 

C = E Pg (s„e- ST (£ J T e *{B u -B T ) +{ r-8 + 4){u-T) du 1) + 

= E(S e- ST (± ^ e „ (Wu -W T ) + (r-S+4)^-T) du „ 1)J 

= Efe-^^T-^o) 



where £t is the process defined by \1J$ but with 7 = r — S + ^- . We see therefore that the problem simplifies 
to the fixed strike Asian pricing problem. 

Let us write down the stochastic differential equation that rules the process (£t)tefo,Tl ■ Using Ito's lemma, 
we get 

f d£ t = + 6 (adW t + (7 + 4)*) 

{ Co = ^o- 

Note that we are faced with a singularity problem near because of the term ^°~^ f . We are going to reduce 
its effect using another change of variables. 
Using Ito's lemma, we show that 

Co = E(e- rT f(S e XT )) (15) 
where Xt = log(£t/£o) solves the following SDE 



dX t = adW t + jdt + 6 l~ x dt 
X = 0. 



-x t _ , 

t~ dt (16) 



Lemma 8 — Existence and strong uniqueness hold for the stochastic differential equation H6\) . 

Proof : Existence is obvious since we have a particular solution X t . The diffusion coefficient being 
constant and the drift coefficient being a decreasing function in the spatial variable, we have also strong 
uniqueness for the SDE (see the proof of Proposition [4]). □ 

Because of the singularity of the term - — j— - in the drift coefficient, the law of (X t )t>o is not absolutely 
continuous with respect to the law of {o~W t )t>o- That is why we now define (Z t ) t > by the following SDE 
with an affine inhomogeneous drift coefficient : 

dZ t = adWt+jdt-^dt . 
Z = X a = 0. 

The drift coefficient exhibits the same behavior as the one in (1161) in the limit t — > in order to ensure the 
desired absolute continuity property. It is affine in the spatial variable so that (Z t )t>o is a Gaussian process 
and as such is easy to simulate recursively. 
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Lemma 9 — The process 

Z t = jJ o *dW s + It (18) 
is the unique solution of the stochastic differential equation J i7| ). 

Proof : Using Ito's Lemma, we easily check that Z t given by (fT5|) is a solution of p7|) . Again, constant 
diffusion coefficient and decreasing drift coefficient ensures strong uniqueness. □ 



Remark 10 — For the computation of the price Co = E (e rT (Soe XT — -?0+) of a standard Asian call 
option, the random, variable e~ rT (Soe ZT — K) + provides a natural control variate. Indeed, since Zt is a 
Gaussian random variable with mean and variance one has 



E (e~ rT (S e z ^ - K)+) = S a e^+^ T M (d + a^T^j - Ke~ rT Af(d) 

where Af is the cumulative standard normal distribution function and d = 2 T _ 

Notice that in Kemna and Vorst \1S^ . the authors suggest the use of the control variate 
rT (^Sq exp ^ oWt + jt dt^j — K^j which has the same law than e~ rT ( K Soe ZT — K) + 



e 



— I crWt + "ftdt is also a Gaussian variable with mean ZT and variance 

T J 2 3 



In order to define a new probability measure under which (Z t )t>o solves the SDE (|16j) . one introduces 



L t = exp 



Z °~l + Z, .„. 1 /•' (> 7 \ ■ 



-,\\\\ --J ( — ) ,1s 



Because of the singularity of the coefficients in the neighborhood of s = 0, one has to check that the integrals 
in L t are well defined. This relies on the following lemma 

Lemma 11 — Let e > 0. In a random neighborhood of s = 0, we have 

\Z S \ < cs2- £ and \X„\ < c S i~ e 
where c is a constant depending on o~,~{ and e. 

Since Ve > 0, 

vz<c^,( e " z ; s 1+z ) 2 <c s - 4 -, 

we can choose e < \ to deduce that L t is well defined. 

,(3t)1? 

Proof: We easily check that the Gaussian process (-Bt)te[o,T] defined by B t = I sdW s is a standard 

Jo 

Brownian motion. Thanks to the law of iterated logarithm for the Brownian motion (see for example Karatzas 
and Shreve [T^] p. 112), there exists ti((d) such thalQ, 

V< < ti(w),|B t (w)| < 



u) is an element of the underlying probability space Q. 



15 



Therefore, 



V* < (3ii(w))4, \Z t {u)\ = IjB^uj) + lt\ < ~^t^ + It. 

' t 3 2 1 32 _ 3 I 



Taking c = max( i CT , Tj) yields 



32-3- 

V* < (3ti. (w)) 3 A 1, |Zi(w)|<cd 



On the other hand, recall that X t = log(&/£ ) = log f ie^*"^* e~ aW ^ u du^ . So, using the law of 
iterated logarithm for the Brownian motion, we deduce that there exists *2(w) such that 

V* < t 2 (u), < - e aWt ^ + ^ [ e- aWu{ul) - lu du < I e <Tti_£+7t / e au ^~ lu du. 



* Jo t Jo 



Denote g(t) = \e at ^ +7 * L e CTtl ' ~ lu du and let us investigate the order in time near zero of this 
function. We have that 

e^'+T* = 1 + a*i" e + 0{t l ~ 2e ) 



t 

I - e 



hence 



<?(*) = 1 + (a + ^— + ©(i 1 -^), 



2 - ( 

so X t (w) < log (<?(*)) ~ (cr + "^r - ~ — )^~ e , which ends the proof for X t . □ 
t-»o I - e 



Proposition 12 — (-^t)ie[o,T] * s a martingale and, consequently, for all g : C([0,T]) — > IR measurable, the 
random variables g((X t )o<t<T) and g{{Z t )o<t<T)LT o,re simultaneously integrable and then 



~-{9{{ x t)o<t<T)} = ^(g{{Zt)o<t<T)L T y 



Proof : The proof is similar to the proof of Proposition 3) 

We have already shown existence and strong uniqueness for both SDE (fT()]) and (fT7]) . Showing that the 
stopping time 

f /"* ( e~ Ys - l + ^\ 2 1 
r n (y) = inf < * e R + such that / ) ds > n \ , with the convention inf{0} = +00, 

[ Jo V °s J J 

have infinite limits when n tends to +00, Qx and almost surely, follows from the previous lemma. 



One has 

Lt = exp 



T e" z < -1 + Z t .„ f T e- z * -1 + Z t (e- z < - 1 + Z t Z t 
aH dZt ~ i aH 2* + 7-y]* 
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Set A(t, z) 



-. The function A : ]0, T] 



is continuously differentiable in time and 



twice continuously differentiable in space. So, we can apply Ito's Lemma on the interval [e, T] for e > : 



A(T,Z T ) = A(e,Z £ ) 



T e -z t 



1 + Z t 



dZ t 



l-Z t + ^--e 



-z% 



-dt 



1 - e~ Zt 
2t 



dt 



Using the lemma |9j we let e — > to obtain 

rT e- Zt -1 + Z t 



A(T, Z T ) 



dZ, 



T 1 _ Zf + ^_ e -z t 



2-I-2 



-dt 



1 - e~ z * 
2t 



dt. 



Then 



Lt = exp 



A(T,Z T )- [ 4(t,Z t 
Jo 



)dt 



where <j) is the mapping 



4>(t,z) 



l + z- 



2+2 



+ 



1 -e 



21 



— I + z / e z — 1 + z z 
( 2t + 7 ~ t 



By (fl5|) and Proposition [12j we get 

C = E[e- rT f(S e ZT )exp 



A(T,Z T )~ [ <f>(t,Z t )dt 
Jo 



(19) 



(20) 



Since for each t > 0, lim <f>(t, z) = +oo and lim </>(t, z) = — oo, it is not possible to apply the exact 

z — ^ — oo z— >-+oo 

algorithm. One can use the unbiased estimator, at least theoretically, if there exists a random variable cz 
measurable with respect to Z such that 



< oo. 



£ (e A{T ' ZT) - [r+Cz)T \f(S e ZT )\eS° Wz-4>{t,z%)\dt^ 
Unfortunately, this reinforced integrability condition is never satisfied : 



Lemma 13 — Assume that f is a non identically zero function. Letpz andqz denote respectively a positive 
probability measure on N and a positive probability density on [0,T]. Let N be distributed according to pz 
and (J7j)ieiM* be a sequence of independent random variables identically distributed according to the density 
qz , both independent conditionally on the process (Zt)t^[o,T]- Then the random variable 



^A(T,Z T )-rT 



1 N 



-(j)(Ui, z Ut 
pz{N)Nlf = \ q z (U z ) 



(21) 



is non integrable. 

Proof : By conditioning on Z, one has 

A ir~ ( e A ^ T - z T)- T \f(S e z T )\ t~tN \<j>(U„Zu t )\ 

Pz(N)N! 1L=1 qz(U t ) 



E (e A ( T ' Z ^- rT \f(S e z ^)\efo m^ t )\dt 



> E ( e A{T,Z T )-rT mSQe Z T) \ e lT W*-*)!* 
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One can easily show that, Vz < and Vt € [-j, T], <p(t, z) > <j>{z) where 



(z) = 



o*(f)» 



— 1 + Z , 2 

y +7+ -2- 



Since 0(z) ~ 2— — -, there exists c < such that for all z < c, <j)(z) > . Hence, 
— oo <j T a 

A > E(' e ^( r ^)-'- T |/(5oe^)| e ^ / l e "^ 1{Zt<o}dt 
> E fe A ( T < z ^- rT \f(Soe ZT )\e-^e^ J % e ~ 2Ztdt 

Using Jensen's inequality we get 

We have seen in the proof of lemma [Til that Z t — %-B t 3 + where (B t )t>o is a standard Brownian motion. 
So, conditionally on Zt, Jt Z t dt is a gaussian random variable and hence A = +oo. 



We are in a situation where e A ^ T ' z ' r '>~ rT \f(Soe ZT )\ E 



N -${Ui,Zu,) 



grable while e A ^ z ^ rT \f(S e ZT ) 

E (e- rT \f(S e z r)\ exp [a(T, Z t ) - </>(t, Z t )dt 
for a given neN*, the random variable 



pz{N)N\ 
pz(N)N\ lli=l q z (Ui) 



nr =1 

(Zt)te[o,T\ 



qz(Ui) 



(Zt)te[o,T] 



is non intc- 



is integrable since 



< oo. Then, a natural idea would consist in considering, 



n ^ — ' 



n 



where (Nj)i<j< n are independent variables having the same law as iV and ( (ff 



(2t)te[o,T] 



are independent 

1 i<j'<n 

sequences having the same law as (Ui) ie ^,, both independent conditionally on the process {Zt)te[o,T]- The 
following general result tells us that this is not sufficient to circumvent intcgrability problems. 



Lemma 14 — Let Y and Z be two real random variables and g : R — > R a given measurable function. 
Assume that g(Z)E(Y\Z) is integrable while g(Z)E(\Y\\Z) is non integrable. Then, when {Yi)\<i< n is 
a sequence of independent random variables having the same law as Y, Vra G N*, the random variable 
g(Z)E (|~ Y2i=i Yi\ \Z) is non integrable. 



Proof : Denote by e, e± and e n three functions satisfying 

VzeR, e(z) = E(Y\Z = z) , e\(z) = E (|Yi| \Z = z) and e n (z) 



1 " 

71 ^ ' 



\Z = z 



On the one hand, since L \g(z)\ \e(z)\Pz(dz) < oo and L \g(z)\ ei(z)P z[dz) = +oo , where Pz is the law of 
Z, we have that J R \g(z) \ ei(z)t {ei ^> 2 \ e ( z )\}P z(dz) = +oo. 
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On the other hand, Vz £ R, 



e„{z) > - 



1 

n 

> I 



{V2<i<n,Kj>0} 



Z = z\+£ 



{V2<i<n,y 3 <0} 



Z = z 



E(F 1 +|Z = ^) P(F a > 0|Z = 



z)"" 1 + E (Yf|Z = z) P (Yi < 0\Z = zf- 1 



ei(z) + e(z) P (Ya > 0|Z = zf _1 + ei(z) 2 " e(z) P (Yi < 0|Z = z)"" 1 



> A 



gi(f) 
4 



l{e l( z)>2| eW |}P (Y > Q\Z = z) n 1 + 2iMl {ei(z) > 2 | e(z) | } P (Y < 0|Z = zf" 1 



> 



l^l{e 1 ( Z )>2| e (z)|} 



Hence, E[g(Z)E(|iEr=i^| 1^)] = / R \g(z)\e n (z)P z (dz) = +oc. 



There is still hope yet. In the proof of Lemma 1131 we saw that integrability problems appear when 
Z t takes large negative values so that (f>(t, Z t ) tends rapidly towards +00. Since lim <f>(t, z) = —00, one 

Z— 5- + 00 

possible issue is to split the function (j>{t, Z t ) into a positive part and a negative part. The first term can be 
handled by the exact simulation technique whereas the second term, which as we shall see in the following 
section presents no integrability problems, can be handled by the unbiased estimator technique. 



2.2.1 An hybrid pseudo-exact method 

We rewrite (|20[) in the following form 



C = E (e A ( T > z ^- rT f(S e z T)efo T *-&z t )4t e -£ 4> + Wt)dtj 



(22) 



Let pz and qz denote respectively a positive probability measure on N and a positive probability density on 
[0,T]. Let N be distributed according to pz and {Ui)i^- be a sequence of independent random variables 
identically distributed according to the density qz, both independent conditionally on the process (Zt)te[a.T]- 
Note that, since e A ^ T ' Z ^~ rT f{S a e z ^)e^ \4>-«,z t )\dt e - 4> + (t,z t )dt = e A{T,z T )-rT f( So e z T) e - So ' *(*,■&«)«« fa 
intcgrable, one has 

C = E\e Ai - T ^- rT f{S e z -) 



N 

n 



PZ {N)N\ \AA q z (Ui) 



di 



(23) 



Remark 15 — There is no hope that this estimator is square integrable. Indeed, one can show as in Lemma 
mUthat E (e^^ ( *' Zt) ) dt ^j = +00 since (4>~(t,z)) 2 is of order z 4 for large positive z. 

The idea then is to apply the exact simulation technique to simulate an event with probability e~ fo < ^ + (*>' z *) d *. 
Since for each t > 0, lim <fi + (t, z) = +00, one needs to bound from above (j> + (t, z), uniformly with respect 

z— > — 00 

to t 6 [0, T], for z > c where c < is a given constant. Thanks to the following lemma, it is possible to do 
so but only uniformly with respect to t € [e, T] for all e > : 



Lemma 16 — For all <t<T, 



sup^,z)<2! + -L + lQ-J_ 
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and 

Ve < 0, ,) < S^±f (, + 7+t , + fe^-il! - -fl. 



Proof : Let > 0. It is useful to distinguish two cases according to the sign of 7 : 
1. 7 > 

We rewrite <j> in the following form 

A e-*-l + z-£ 1-e-' (I 7 V 7* ^-(^l) 2 (e~» - l) 2 - (z A l) 2 



cr 2 t 2 t V 2 o- 2 J o-H 2a 2 t 2 2a 2 t 2 

First note that e ^~J+r^ < 0, ±=f^ (I - ^) < I (I - 4,) + and ^""^i** 1 * 2 < 0. Moreover, 
70 z 2 -(zAl) 2 1 / z 1 /z\ 2 (f A±) 2 



er 2 t 2o~ 2 t 2 a 2 y't 2 V i 

< I f { 

5t otherwise 



Consequently, z) < £ + ^ + | (§ - ^) + . 

2. 7 < 

Now we rewrite </> in the following form 



e-*-l + z-4r e^-l + z ( e - z - l) 2 - z 2 . 1 - e 
0(*> *) = 575 + 7" 



<7 2 i 2 ' a 2 t 2a 2 t 2 2t 

It is then easy to show that <j) + {t, z) < ^. 

Note that ^ < ^y + ^ + j (| — ^) + ■ Hence, gathering the two cases yields the first part of the lemma. 
Let now z e [c, 0] for a given negative constant c. We rewrite </> in the following form 

, e-'-l + z, + , (e- z -l) 2 z 2 l-e- z e~ z -l + z 



>0 for z<0 <0 for z<0 



Since 9 Z 



e~ z -l + z^ , , (e- z -l) 2 z 
(l + 7 + f) + - s 



<7 2 t 2 ' y 2cr 2 i 2 er 2 i 2 

z < 0, one has that 



z ^ 2 2 ' 1- e- 2z -2z + t~f+(l- e- z ) . . r 

is negative for all 



i 2 <7 



2 



sup </>+(t,z) < 5-5 (l + 7 +t)+ V 



,e[c,o] <^ 2 ' ' 2a 2 t 2 <7 2 t 2 ' 



□ 



This lemma suggests to apply the exact algorithm on [e,T] for a fixed positive threshold e. It remains 
to handle the time interval [0,e[. Thanks to the following lemma, we that (j> + (t,Z t ) can be approximately 
bounded from above for small t, almost surely, by a function of t. The idea is then to extend the exact 
simulation algorithm by simulating an inhomogeneous Poisson process. Of course, this hybrid method is no 
longer exact since the positive threshold for which the upper bound holds is random. 
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Lemma 17 — For all rj > 0, there exists a random neighborhood oft — O such that 



o-(/.z,)r : (fj + f )' " (.-'4) 



Proof : We rewrite (1191) this way 

,, N /l-e" z e" z -l + z\ 1 (l-z + ^-e- z -Ue~ z -l + z){e- z -l-z)\ 1 



2 ' a 2 J t \ a 2 J t 2 

and make the following Taylor expansions 

l-*+$-e-'-h(e-'-l + z)ie--l-z) = 2_ + 
a 1 3a 2 

, l-e~ z e~ z -1 + z 1 2 , 
and + 7 - = -z + 0(z 2 ). 

On the other hand, we have seen in the proof of lemma [TT] that there exists a random neighborhood of zero 
such that Z t < ct2~ v where c = max( i a _ R , We conclude that, in a random neighborhood of zero, 

37 3 

□ 



2.2.2 Numerical computation 

For numerical computation, we are going to use the following set of parameters : Sq = 100, K = 100, 
a = 0.2, r = 0.1, 5 = and T = 1. To fix the ideas, let us consider a call option. The price Co writes as 
follows 

Co = E (^MT.Z T )-rT _ K f ^ Cp rO^M j e -/ o T 0+(^ t )^ . 

where iV ~ 'P(cp) and (J7i)i>i is an independent sequence of independent random variables uniformly dis- 
tributed in [0,T]. The parameter c p > is set to one in the following. We give a description of the hybrid 
method we implement : 

Algorithm 2 

On the time interval Ij := [t^ttj 

1. Simulate Z t ,Zt and a lower bound mi for the minimum of (Zt)tei- (use the fact that Zt = 
fB t 3 + %t where (B t )t>o is a standard Brownian motion). 

2. Find M j > such that Vi G Ij,(j) + {t, Z t ) < M j (use LemmaUb}). 

3. Simulate an homogeneous spatial Poisson process on the rectangle Ij x [0, AP] and accept (respectively 
reject) the trajectory simulated if the number of points falling below the graph {<j) + (t, Zt))t£ij is equal 
to (respectively different from) zero. 
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Carry on this acceptance rejection algorithm until reaching a time interval Ij for a chosen J€ N*. On the 
remaining time interval [0, 2 i+i L use the same acceptance/rejection algorithm but with an inhomogeneous 
spatial Poisson process this time (use Lemma \T7^ . 

In table[3j we give the price obtained by our method for different values of the positive threshold e = ^prr- 
The number M of Monte Carlo simulations is equal to 10 5 and the true price is equal to 7.042 (computed 
using a Monte Carlo method with a trapezoidal scheme and a Kemna-Vorst control variate technique). 





Price 


CPU 




6.9394 


7s 


f - T 

t — 9T 


6.9590 


10s 


e _i 

t — 2 6 


6.9703 


13s 


f - T 

6 — oF 


6.9952 


17s 


e - 5TTT 


7.0423 


21s 



Table 3: Price of the Asian call using the hybrid-pseudo exact method. 

Clearly, the method is not yet competitive regarding computation time. Nevertheless, unlike the usual 
discretization methods, it is not prone to discretization errors. 

3 Conclusion 

In this article, we have applied two original Monte Carlo methods for pricing Asian like options which 
have the following pay-off : (aSx + /3 S t dt — K) + . In the case a ^ 0, wc applied both the algorithm of 
Beskos et al. [T] and a method based on the unbiased estimator of Wagner [33] and more recently the Poisson 
estimator of Beskos et al. [2] and the generalized Poisson estimator of Fearnhead et al. [Bj . The numerical 
results show that the latter performs the best. The more interesting case a — 0, which corresponds to usual 
continuously monitored Asian options, can not be treated using neither the exact algorithm, nor the method 
of exact computation of expectation but we investigate an hybrid pseudo-exact method which combines the 
two techniques. More generally, this hybrid method is an extension of the two exact methods and can be 
applied in other situations. 

From a practical point of view, the main contribution of these techniques is to allow Monte Carlo pricing 
without resorting to discretization schemes. Hence, we are no longer prone to the discretization bias that 
we encounter in standard Monte Carlo methods for pricing Asian like options. Even though these exact 
methods are time consuming, they provide a good and reliable benchmark. 
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4 Appendix 



4.1 The practical choice of p and q in the U.E method 

The best choice for the probability law p of N and the common density q of the variables (Vi)i>i 
is obviously the one for which the variance of the simulation is minimum. In a very general setting, it is 
difficult to tackle this issue. In order to have a first idea, we are going to restrict ourselves to the computation 



U(JV)JV!ll«(K) 



where g : [0,T] -> DR. 



Lemma 18 — When g is a measurable function on [0, T] such that < / \g(t)\dt < +oo, the variance of 

Jo 

1 N 
vTa/T IT 



g(Vi) . 



p{N)N iii q{Vi) 
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(lo\9(t)\ 



dt 



exp 
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Proof : Minimizing the variance in ([7]) comes down to minimizing the expectation of the square of 
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Using Cauchy-Schwartz inequality we obtain a lower bound for F(p, q) 
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n=0 
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= <*P^2jf \9{t)\dt 
We easily check that this lower bound is attained for q opt and p op t- 



The optimal probability distribution p opt is the Poisson law with parameter / \g(t)\dt. This justifies our 

Jo 

use of a Poisson distribution for p. 
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4.2 Simulation from the distribution h given by ( 1131) 



Recall that 



h(u) = Cexp A(u) 



(u - X f 



2T 



6 exp — u H — (1 

(7 CT 



2T 



where C is a normalizing constant. 
The expansion of the exponential e~ 



at the first order yields 



, , . r ,7 ^ PS (u - X f 
h{u) ps (7 exp ( —u H u — 

1 (7 (T 2i 



= C exp — 



2T 



This suggests to do rejection sampling using the normal distribution with mean Xq + t (t+ /3s °) an d 
variance T as prior. Unfortunately, for a standard set of parameters, this method gives bad results. Even a 
second order expansion of e~ au which also modifies the variance does not work. 

In order to get round this problem, we evaluate the mode u* of h. We have 



1 a a 



So, h'iu*) = if and only if 



T 



1 , PS _ c 



7 * 
exp [ — u 
a 



a 2 [ 



) 



(U* ~ Xq) 2 

2T 



U* - Xq 

T 



= 



which writes 



a{u* -X - lT)e a ( u '- x °-« T ) = TpS e- aX °^ T . 



The function x i— > xe x is continuous and increasing on [0, +oo[ and so is its inverse which we denote by 
W . Since T f3Soe~ aX °~ lT > 0, we deduce that h is unimodal and that its mode satisfies 

7 T + W (/3S Te-~< T -' 7Xo ) + aX 

u = . 

a 

The function W is the well-known Lambert function, also called the Omega function. It is uniquely valued 
on [0, +oo[ and there are robust and fast numerical methods based on series expansion for approximating 
this function (see for example Corless et al. [4]). 

Numerical tests showed that performing rejection sampling using a Gaussian distribution with variance 
T and mean u* instead of Xq + T (' 7 + /3So ) gives plain satisfaction. In table |H we see that for arbitrary choice 
of the parameter ^-g, the acceptance rate of the algorithm is always high (of order 70%) and that the 
computation time is low. 



a 


Nb of simulations 


Acceptance rate 


Computation time 


0.2 




61% 


3s 


0.5 


10 6 


68% 


3s 


0.8 




80% 


2s 



Table 4: Acceptance rate of the rejection algorithm of simulating from the distribution h in (|13|) with 
S a = 100, a = 0.3, T = 2 and r = 0.1. 
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